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Abstract. We study the gas accretion onto a supermassive black hole (SMBH) using 
the 3D SPH code GADGET-3 on scales of 0.1 - 200 pc. First we test our code with 
spherically symmetric, adiabatic Bondi accretion problem. We find that our simulation 
can reproduce the expected Bondi accretion flow very well for a limited amount of time 
until the effect of outer boundary starts to be visible. We also find artificial heating of 
gas near the inner accretion boundary due to the artificial viscosity of SPH. Second, 
we implement radiative cooling and heating due to X-rays, and examine the impact of 
thermal feedback by the central X-ray source. The accretion flow roughly follows the 
Bondi solution for low central X-ray luminosities, however, the flow starts to exhibit 
non-spherical fragmentation due to thermal instability for a certain range of central L x , 
and a strong overall outflow develops for greater L x . The cold gas develops filamentary 
structures that fall into the central SMBH, whereas the hot gas tries to escape through 
the channels in-between the cold filaments. Such fragmentation of accreting gas can 
assist in the formation of clouds around AGN, induce star-formation, and contribute to 
the observed variability of narrow-line regions. 



1. Introduction 

Feedback from active galactic nuclei (AGN) is considered to be one of the important 
processes in regulating galaxy growth over cosmic time. Several research gr oups have 



implemented AGN feedback in cosmological hydrodynamic simulations (e.g. lDi Matteo et al 



20051 : iBooth & Sc haye 2009), however they all have to assume a subgrid model for the 



accretion onto SMBH and feedback with uncertain parameters due to limited resolution 
in cosmological simulations. 

Here we take a different approach and perform small-scale simulations of gas ac- 
cretio n onto SMBH using the GADGET-3 SPH code (originally described by lSpringell 



20051) on small-scales of r < 200 pc. Our long-term goal is to develop a better model of 
AGN feedback for large-scale cosmological simulations based on the results of small- 
scale BH accretion simulations. In this article, we report our initial results of adiabatic 
Bondi accretion simulations, a s well as those with radiative cooling and heating by X- 
rays from the central SMBH (iBarai et alJ 1201 lbh . We run a large set of simulations 



of spherically symmetric Bondi accretion with different resolution, volume, and initial 
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Figure 1 . Radial profile snapshot of gas particles in a run with X-ray heating and 
cooling, showing radial velocity, density, temperature and Mach number. The solid 
line is the median value, the dashed lines are the 95 percentile, and the dotted lines 
are the min/max values. The red curve in each panel shows the free-fall scaling, and 
the blue curve shows the corresponding ZEUS simulation results: vzeus and Vff are 
indistinguishable, Pzeus is at a constant offset from p$ because of different normal- 
ization. The green curve in the bottom-left pane l indicates the free-f all temperature 
with only adiabatic processes. Figure taken from lBarai et al.l (l2011bl) . 



conditions. The Bondi problem is ideal for testing any hydrodynamic code, since it has 
a semi-analytic solution. 



2. Simulations of Bondi Accretion onto a SMBH 

The central SMBH has a mass of Mbh = 10 8 M o , and the initial conditions are gener- 
ated using y = 1.01, poo = 10" 19 g/cm 3 , and T^t = = 10 7 K. The corresponding 
Bondi radius is Rb = GMbu/c^ ^ = 3.0 pc, the Bondi time is ts = Rb/c s = 7.9xl0 3 yrs, 
and the sonic radius is R s = 1.5 pc. The total initial mass of the gas sphere is in the 
range of 10 5 - 10 7 M o depending on the setup, and the particle counts were varied be- 
tween 64 3 - 256 3 for resolution test. The inner boundary radius is ri n = 0. 1 pc from 
the central SMBH, and we regard that the particle is accreted to SMBH once it crosses 
ri n . The outer boundary radius was varied between r ni1t = 5 - 200 pc. For the full list of 
runs and parameters, see Table 1 of iBarai et al. I (l2011bh . 



First we find that the adiabatic simulations reproduce the expected Bondi solution 
very well for a limited time. Since the gas particles near the outer boundary escape due 
to adiabatic expansion, after a certain time the mass accretion rate start to decrease from 
the Bondi rate. This is a problem owing to the difficulty of setting a boundary condition 
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in the SPH method, and it does not arise with the mesh codes, as one can simply fix the 
outer boundary with and Too- 

In the second set of runs with rad i ative c ooli ng and heating by X-rays, we follow 
the approximate treatment of iBlondinl (1 19941 ) and IProgal (120071 ). assuming that the 10 
keV bremsstrahlung radiation from the central SMBH is illuminating the optically thin 
gas. In Figure [H we show the particle radial profile for various quantities in a run with 
L x = O.OlLEdd* where LEdd is the Eddington luminosity. We find that the simulation 
roughly follows the Bondi flow solution even with the X-ray heating and cooling, except 
at r < 0.5 pc where we see a spurious heating due to artificial viscosity (AV) of SPH 
particles. We have confirmed that this effect is due to AV by turning it on and off. With 
no AV, the overheating does not occur, and the gas particles overshoot the SMBH due 
to lack of viscosity and many are not accreted. The mass inflow rate at r m in this run is 
enhanced over the Bondi accretion rate by a factor of a few due to cooling. 



3. Non-spherical Fragmentation due to Thermal Instability, and Outflows due to 
X-ray Thermal Feedback 

To study the effect of X-ray thermal feedback, we gradually increase Lx/L^dd from 
5 x 10" 5 to 5 x 10" 2 . In these runs, the initial condition contains 12.7 million particles 
for a gas sphere of 9.77 x 10 6 M o , each gas particle mass of 0.791M o , adiabatic index 
y - 5/3, r out = 200 pc, = 10" 23 g/cm 3 , and = 10 5 K. The Bondi radius for these 
parameters is 183.9 pc, therefore most of our computational volume is within the Bondi 
radius and we fully resolve below the Bondi radius with minimum smoothing length 
of ~ 0.1 pc. In general we find that, with increasing Lx, the mass accretion rate at r m 
decreases, and the outflow rate at r out increases. 

When L x becomes > O.OlLEdd? we start to see an interesting transition from inflow 
to outflow, and non-spherical fragmentation of gas into multiphase medium takes place 
due to thermal instability. The filamentary cold gas continues to flow in, and the hot 
gas tries to escape through the channels between cold filaments. Examination of the 
flow motion shows that the filaments get stretched, fragments, and the 'clouds' merge. 
Such fragmentation of accreting gas can assist in the formation of clouds around AGN, 
induce star-formation, and contribute to the observed variability of narrow-line regions. 

Figure [2| shows an example of such multiphase structure in a run with L x = 
O.OlLEdd? when the outflow is still not very prominent and the hot gas is still being 
accreted to SMBH. Studies with ID & 2D ZEUS simulation shows that this fragmen- 
tation does not occur in a spherically symmetric Bondi flow, unless some perturbation 
is introduced by hand. The SPH simulations inherently have tiny fluctuations in the 
density field due to its algorithmic nature, which can be amplified through thermal in- 
stability. 

As L x is increased to >0.01LEdd> the outflow starts to dominate over the inflow, 
and the hot gas escapes through the channels between cold filaments, as shown in the 
left panel of Figure [3] We find that the transition from inflow to outflow occurs in- 
between Lx/L^dd = 0.01 - 0.02, but note that this transition luminosity would depend 
on the value of p^. In other words, the more relevant parameter for the transition is 
the range of photoionization parameter £ oc L x /p for the unstable branch of the T - £ 
equilibrium curve. 

With a high enough L x (= 0.05LEddX a strong outflow is produced due to strong 
thermal feedback, and the mass outflow rate at r - r out increases dramatically. A jet- 
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Figure 2. Density and temperature cross-sections in a run with Lx/^Edd = 0.01. 
The top two panels show the inner ±30 pc, and the bottom two the inner ±4 pc of the 
[x - y] plane through z = 0. Cold, dense filamentary structure has developed due to 
thermal instability, which is falling into the SMBH rapidly. The hot gas is trying to 
escape, but a strong outflow has not developed yet due to low L x . 



like plume of hot and buoyant gas escapes from the central region to the outer boundary, 
as shown in the right panel of Figure [3] 

T hese results on the non-spherical flow have been reported in our second paper in 
series dBarai et al.ll2011ah . where we examined the T - £ equilibrium curve in detail and 
determine the unstable ^-range. 



4. Conclusions 

We find that GADGET-3 SPH code can reproduce the spherically symmetric Bondi 
accretion flow properly with mainly two limitations: 1) the gas particles escape from 
the outer boundary due to adiabatic expansion; 2) spurious heating is observed near the 
inner boundary around SMBH due to the artificial viscosity of SPH. 

We also examined the impact of radiative cooling and heating due to X-rays. The 
accretion flow roughly follows the Bondi solution when the central X-ray luminosity 
is relatively low, but the mass accretion rate at r m is enhanced by a factor of a few 
over the Bondi accretion rate due to cooling. The simulation starts to exhibit non- 
spherical fragmentation due to thermal instability once L x exceeds ^ O.OlLEdd- At 
L x = 0.02LEdd> the hot gas escapes in-between the cold filaments. When L x is further 
increased to 0.05LEdd> the outflow completely dominates over the inflow, and most 
of the gas escapes from the computational volume. A jet-like outflow feature is also 
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Figure 3. Temperature cross-section in runs with Lx/L^dd = 0.02 (left) and 0.05 
(right). The left panel shows the inner ±100pc, and the right panel shows the entire 
computational volume of ±200 pc. In the left panel, the cold, dense, filamentary 
structure has developed due to thermal instability, which is falling into the SMBH. 
The hot gas tries to escape through the channels in-between the cold filaments. In 
the right panel, the outflow dominates over the entire computational volume, and the 
hot plume of gas is escaping from the center to the outer boundary. 

observed at this L x , which selects a preferential direction for hot gas to escape rapidly. 
Such non-spherical features of accreting gas can assist in the formation of clouds around 
AGN, induce star-formation, and contribute to the observed variability of narrow-line 
regions. 

In the future, we plan to implement rotation, radiation pressure, and different initial 
geometry of gas distribution. We will measure the efficienci es of thermal, radiative , and 
kinetic feedback, and compare them with those measured by lKurosawa et al.1 (l2009h and 
those used in the cosmological simulations. 
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